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This paper extends our earlier studies of free energy functions of density and crystalline order parameters 
for models of supercooled water, which allows us to examine the possibility of two distinct metastable liquid 
phases [J. Chem. Phys. 135, 134503 (2011) and arXiv: 1107.0337 1^2]. Low-temperature reversible free energy 
surfaces of several different atomistic models are computed: mW water, TIP4P/2005 water, SW silicon and 
ST2 water, the last of these comparing three different treatments of long-ranged forces. In each case, we show 
that there is one stable or metastable liquid phase, and there is an ice-like crystal phase. The time scales 
for crystallization in these systems far exceed those of structural relaxation in the supercooled metastable 
liquid. We show how this wide separation in time scales produces an illusion of a low-temperature liquid-liquid 
transition. The phenomenon suggesting metastability of two distinct liquid phases is actually coarsening of 
the ordered ice-like phase, which we elucidate using both analytical theory and computer simulation. For 
the latter, we describe robust methods for computing reversible free energy surfaces, and we consider effects 
of electrostatic boundary conditions. We show that sensible alterations of models and boundary conditions 
produce no qualitative changes in low-temperature phase behaviors of these systems, only marginal changes 
in equations of state. On the other hand, we show that altering sampling time scales can produce large and 
qualitative non-equilibrium effects. Recent reports of evidence of a liquid-liquid critical point in computer 
simulations of supercooled water are considered in this light. 
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I. INTRODUCTION 

This is our second paper examining whether molecular 
simulation provides support for the hypothesis that su- 
percooled water possesses two distinct liquid phases with 
a reversible coexistence line ending at a critical point. ^ 
In the first (Paper I) ^ we described our results for perti- 
nent free energy functions of three different models: the 
mW model of water,l21 a variant of the ST2 model for 
water,S5 and the Stillinger- Weber model for Each of 
these models has an equilibrium liquid phase with pair 
distributions and thermodynamic anomalies like those of 
water, and each has an equilibrium phase transition like 
that of water-ice freezing. For each, others have claimed 
numerical eviden ce o f liquid-liquid transitions at super- 
cooled conditions but our calculations described in 
Paper I reveal no such behavior. Rather, for each sys- 
tem we found that there can be at most one stable or 
metastable liquid phase, and this liquid can coexist with 
crystalline ice. Here, we establish that results contrary 
to our findings derive from lack of equilibration of slow 
fluctuations in long range order. We also present new 
calculations for several additional models reaching consis- 
tent conclusions in each case. Specifically, for computer- 
simulation models of water that exhibit liquid and ice-like 
phases, there is no liquid- liquid transition, but there is 
non-equilibrium coarsening of ice that others have mis- 
interpreted as evidence of a liquid-liquid transition. 
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Figure [T] shows the relevant part of water's phase dia- 
gram and corresponding free energy surfaces. The liquid 
is the stable equilibrium phase for temperatures above 
the melting temperature, i.e., T > T^, and it is unstable 
below a stability temperature, i.e., T < Tg. In between, 
for Ts < T < Tin, the liquid is metastable with respect 
to crystal ice. Throughout much of that intermediate 
region, structural reorganization of water is slow, and it 
becomes slower in a super- Arrhenius fashion as temper- 
ature is lowered.l^ This sluggishness can present prob- 
lems for straightforward molecular simulation, as noted 
below, but it is not so sluggish to prevent certain crys- 
tallization of water when the liquid is cooled close to or 
below Ts- Coarsening of water in that regime occurs on 
the time scale of microseconds - fast for experiment, but 
slow for simulation.!^ All speculations on the existence of 
a liquid-liquid phase transition in water locate that tran- 
sition near or below Tg, the so-called "no-man's land" for 
liquid water. As such, it is difficult for experiments to 
prove or disprove the liquid-liquid hypothesis. It is left 
to simulation, which can reversibly control crystalliza- 
tion, to see if such an idea could be correct within the 
purview of statistical mechanics for plausible models of 
water or water-like systems. 

Calculations of free energy functions of relevant order 
parameters are required when using simulation to estab- 
lish phase behavior.'^S As noted, such calculations can be 
difhcult, especially for supercooled water because fluctu- 
ations in this regime are collective and slow. To address 
this difficulty and sort out the phase behavior of super- 
cooled water, we have found it convenient to consider two 
order parameters. One is molecular density, p, that dis- 
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FIG. 1. Phase diagram, free energy surfaces and typical con- 
figurations of cold water. Typical of all systems considered in 
this paper, the specific pictures render results from molecu- 
lar simulations of one particular model. Quantitative scales 
of temperature, pressure and free energy depend upon simu- 
lation model, and these scales are omitted here because this 
figure serves a qualitative purpose only. For experimental wa- 
ter, the phase diagram covers pressures p ranging from 1 bar 
to 3 kbar, and temperatures T ranging from 150 K to 300 K. 
Tm and Ts stand for temperatures at the melting and liquid- 
stability lines, respectively. Blue region in the p-T plane is 
where liquid is stable, red region is where liquid is metastable 
with respect to ice I, and grey region is where liquid is unstable 
(i.e., it is the liquid's "no-man's land"). Corresponding free 
energy surfaces are shown above as functions of density p and 
crystal-order variable Qe- Corresponding molecular configu- 
rations shown below are cuts through a simulation box at the 
ends of trajectories that are initiated in a liquid configuration 
and that run for times shorter than required to crystalize the 
entire sample. Molecules located in crystal-like regions are 
colored red. 



tinguishes different amorphous phases. The other can be 
the Steinhardt-Nelson-Ronchetti Qq that distinguishes 
an amorphous phase from a symmetry-broken crystalline 
phase.'23[2ll The two variables fluctuate on significantly 
different time scales. For example, the liquid structural 
relaxation time (i.e., the equilibration time for p) around 
T « Ts is of order 10~^ s or shorter, whereas the relevant 
equilibration time for Qq in this regime is the time to 
form a crystal, 10^^ s or longer. This wide separation of 
time scales is typical of systems undergoing crystalliza- 
tion transitions:— In the case of water, we will see in this 
paper that it is a principal source of confusion in simula- 
tion studies that claim evidence for a liquid-liquid phase 
transition. 

To foreshadow this point, consider the equilibrium 
joint distribution function for the order parameters, 
P{PjQ&)- It is related to the free energy (or reversible 
work) surface for these variables in the usual way: 
F{p,Qq) = —AjbT In P(p, Qg), where is Boltzmann's 



constant. This is the free energy function illustrated in 
Fig. [ij Over time scales large compared to those of liq- 
uid relaxation but possibly not large compared to those 
of crystal formation, the joint distribution is in general a 
non-equilibrium distribution, 

Pncip,Q6,t) ^ P{p\Qe) PnciQ6,t) , (1) 

where P{p\Qq) is the equilibrium distribution for p 
given a specific value for Qq, and PnciQe^t) is the non- 
equilibrium distribution for Qq. The non-equilibrium 
distribution depends upon the protocol with which the 
system is prepared, and its time dependence is irre- 
versible. For large enough t, presuming ergodicity, 
PneiQe^t) approaches the equilibrium P{Qq). But this 
limit can require simulation times thousands of times 
longer than those needed to equilibrate p. Not account- 
ing for this behavior can give the illusion of a reversible 
polyamorphism because the non-equilibrium free energy, 
—kBTln[P{p\QQ) PnciQe, t)], can have a \ow-Qq basin for 
times shorter than those required for Qq to diffuse to- 
wards its equilibrium crystal value at high Qq. 

This possibility, which we refer to as "artificial 
polyamorphism," can be appreciated by comparing the 
free energy surfaces shown in Fig.[l] In particular, imag- 
ine studying the system on time scales where Qq can 
diffuse over no more than the left halves of the pictured 
free-energy panels. -PnolQej^) would then be peaked at 
a low value of Qq, even for cases where a high Qq value 
would be the correct equilibrium value. Thus, if Qq is 
limited in this way to small values, the low-temperature 
(i.e., left-most) panel would then yield a pseudo free en- 
ergy, — fceTln Pne(Pj Qe, 0; with an illusory "amorphous 
basin" at a density lower than that of the metastable 
liquid. This irreversible behavior was discussed in the 
Supplement to Paper I.^^ Specifically, we showed that for 
liquid water at pressures and temperatures in or close to 
"no man's land," small values of Qq will survive while 
the crystal phase begins to coarsen. The bottom left of 
Fig. [l] shows a configuration of water in that regime. 

Section II of the current paper provides a quantitative 
theoretical analysis of this behavior. It shows specifically 
how the polyamorphism of Refs. [TOl El [HI HU and [121 is 
an irreversible effect reflecting the time-scale separation 
between fluctuations in p and fluctuations in Qq. Dur- 
ing the time of coarsening, the faster order parameter, p, 
fluctuates between typical crystal and liquid values. Ref- 
erences fT^ and lTCl report this type of behavior, which they 
call "phase flipping." On the time scale of the flipping, 
the drift in Pnc{QQ,t) can be almost imperceptible, but 
drift it does. The authors of Refs. [15| and ^16) describe the 
flipping as evidence of a second metastable liquid. The 
analysis of Section II shows that this flipping is not a 
consequence of such metastability, but rather the coars- 
ening of the crystal from the unstable or nearly unstable 
liquid, occurring steadily and irreversibly on a time scale 
long compared to those considered in Refs. [151 and [T6l 
This finding could already be anticipated from the Sup- 
plement to Paper and from Moore and Molinero's 
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previous study of crystallization of mW waterP^I 

The specific pathways by which simulated models 
coarsen depend upon system size. For example, free 
energy barriers separating coexisting equilibrium basins 
scale as N^^^, manifesting the presence of an interface.^® 
Similarly, with increasing N, the thermally accessible 
widths of the two basins decrease as N~^^'^, and the 
width of the barrier growsf^ Due to the growing barrier 
height and barrier width, the frequency of spontaneous 
events carrying a system between stable phases is there- 
fore vanishingly small in the limit of large N. Similarly, 
in or near the region of liquid instability (i.e., T close to 
or lower than the crossover temperature Tg), the slope 
towards the crystal basin will increase in magnitude as 
N increases. This behavior will affect the rate of "phase 
flipping" during coarsening, giving the transient impres- 
sion of a non-equilibrium barrier between low-density and 
high-density states that changes with TV. 

These finite-size effects are fundamental to the nature 
of phase transitions. Establishing the existence of a phase 
transition requires studying system-size dependence, for 
example, by computing changes in free energy barriers 
with respect to changing N. No such computations have 
yet been performed for putative liquid-liquid transitions 
in models of water that exhibit water-like structure of the 
liquid and crystal phases. To do so requires algorithms 
that can attend to the collective nature of systems under- 
going phase transitions. Free energy methods are among 
the tools that are suitable for the task, provided they are 
combined with trajectory algorithms that are appropri- 
ately efficient and reversible .'23 

In Section III, we detail how pertinent free energies 
can be computed for supercooled water, and we con- 
sider different variants of the ST2 model as applications. 
Juxtaposition of free energy surfaces for three different 
variants indicates that reasonable changes in electrostatic 
boundary conditions do not change general phase behav- 
iors. Section IV presents free energy surfaces obtained 
for other systems: the mW of water, the TIP4P/2005 
model of waterjI^Hland the SW model of Si. In every case, 
the models are found to exhibit one stable or metastable 
liquid phase plus ice-like crystal phases. Coexistence be- 
tween two distinct liquid phases does not occur. A sum- 
mary of our findings is given in Section IV, and Appen- 
dices A, B and C present further details and results. 

Reversibility is particularly important to the issues ad- 
dressed here and in Paper I. Distinct reversible phases 
can be interconverted, with properties that are inde- 
pendent of the paths by which they are prepared. Re- 
versible liquid phases are thus not the same as amorphous 
solids or glasses. The former are reversible and ergodic, 
so their measured stationary behaviors are independent 
of history. The latter, like high-density or low-density 
amorphous ices (HDA and LDA), are not ergodic, so 
their behaviors depend much on history (i.e., prepara- 
tion protocols). Observed transitions between HDA and 
LDA phases, ■^^ therefore, are necessarily different than 
reversible liquid-liquid transitions. Melting amorphous 



TABLE I. Separation of timescales for fluctuations in density 
and long-ranged order. 



Model 



mW 
ST2 
ST2 

Experiment 



10^ MD^ I(FMD^ 



10^ MC^ 
10'' MC^ 
10^ Pd 



10" MC^ 
10** MC^ 
> 10" pJ^ 



^ Liquid structural relaxation time at temperatures close to Tg. 
^ Auto-corrclation time for Qe fluctuations in the liquid at 
temperatures close to Tg. 

NPT molecular dynamics steps at T = 200 K, p = 1 bar.l^l 
d NPT hybrid Monte Carlo steps at T = 235 K, p = 2.2 kbar. 

This work; see text for details. 
" NPT single particle Monte Carlo steps at T = 229 K, p = 2.2 

kbarESl 

^ Estimated for T = 220K and p = 1 bar from analysis in Ref. 1391 
s Estimated from the critical cooling rate of lO^K/s needed to 
form amorphous ice.l^ 



ice to produce a non-equilibrium liquid that then crys- 
talizes is different too^ 

Crystallization following the melting of glas^^H g^^^^ 
crystallization following the rapid quench of water into 
the liquid's "no-man's land"'2S are much like non- 
equilibrium dynamics evolving from low to high Qg on 
the middle and left free energy surfaces pictured in Fig.[l] 
an observation worthy of future study. But these inter- 
esting non-equilibrium processes and the transitions be- 
tween different amorphous solids of water are not our fo- 
cus in this work. Rather, we are concerned with whether 
water-like systems when constrained to not freeze can ex- 
hibit two distinct liquid phases. If such reversible poly- 
morphism were possible, these systems could also exhibit 
a second critical poin t as Stanley and many of his co- 
workers have proposed.Sla[ ^° ' " ' ^'^ ' ^'^ ' ^^ ' 35 3l3Hlif^ instead, re- 
versible molecular simulation models exhibit only ice and 
one liquid, then the symmetry differences between ice 
and liquid exclude the possibility of an associated criti- 
cal point. We believe the systematic evidence provided 
herein and in Paper I indicates that there is only one 
liquid and no low- temperature critical point. 



II. THEORY OF COARSENING AND ARTIFICIAL 
POLYAMORPHISM IN COMPUTER SIMULATION OF 
WATER 

This section provides a quantitative theoretical analy- 
sis showing the difficulty in obtaining correct reversible 
free energy surfaces of supercooled water. We do so by 
examining the effects of time-scale separation for dynam- 
ics on a reversible free energy surface. The particular 
surface we employ is the free energy F(p, Qq) derived in 
Paper I for a variant of the ST2 model. This free energy 
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FIG. 2. Slow relaxation behavior and its consequences for free energy calculations, a) The reversible free energy surface for 
216 molecules with the ST2a potential energy function at temperature T — 235 K and pressure p — 2.2 kbar. (See text for 
definition of the ST2a potential.) Contour lines are separated by 1.5 ksT, and statistical uncertainties are about 1 ksT. b) 
Negative logarithm of the non-equilibrium distribution for crystal order, Qe, as it relaxes from the liquid state. It is computed 
from the Fokker-Planck equation with the free energy surface given in Panel (a) under the assumption that the density, p, 
remains at equilibrium with the instantaneous value of Qg. (c) and (d) Non-equilibrium pseudo free energy surfaces computed 
from Eq. [l]at two intermediate stages of relaxation, t = 10 rq^ and t — 1, 000 tq^. The unit of time, rq^, is the autocorrelation 
time for Qq fluctuations in the liquid basin (i.e., at small Qe)- The reduced density is p = {p — p^ti) / ^p, where pxti is the mean 
density of the crystal basin (i.e., at large Qe), and Ap is the difference between the mean densities of the liquid and crystal 
basins. Contour lines are separated by 1 UbT and statistical uncertainties are about 1 k^T 



is shown in Fig. [2]ja). The methods used to obtain that 
surface are the subject of the next section, but here we 
only need to assume that there is such a surface, and that 
it is qualitatively hke the surface shown in Fig. |2j[a) . 

Free energy surfaces for several models are derived in 
Sections III and IV. The generic features of the free en- 
ergy surfaces are the same for all the models. There is a 
liquid basin at small Qg s,nd large p, and a crystal basin 
at large Qq and small p. For a given temperature, T, the 
relatively stabilities of the basins are controlled by the 
pressure. The free energy at pressure p is related to the 
free energy at pressure p + Ap in the usual way. 



F{p,Qe;p + Ap,T) ^ F{p,Qe;p,T) + ApN/p. (2) 

Accordingly, lowering pressure tips the surface towards 
the crystal basin, and raising pressure tips it towards the 
liquid basin. 

In cases where the crystal is stable but the system is 
prepared in the liquid, an irreversible drift towards the 
crystal will occur. To the extent that p and Qe are the 
principal slow variables, this coarsening can be described 
in terms of motion on the F{p, Qg) surface. By using 
this perspective, and specifically by adopting the free en- 
ergy surface pictured in Fig. , we illustrate here the 
generic behavior of early-stage coarsening of ice. The 
behavior is not specific to the particular free energy sur- 
face. Rather, it is general consequence of a separation 
of time sales, where the density p equilibrates on time 
scales that are at least two orders of magnitude shorter 
than the time scales on which Qq fluctuates. Such sepa- 
rations of time scales are typical in natural and computer 
simulated supercooled water. See Table [ij 



A. Time dependence of Pno(Q6, i) 

Due to the separation in time scales, relaxation of 
PnciQe, t) can be accurately estimated by assuming den- 
sity p is always in equilibrium with the current value of 
Qq . An appropriate Fokker-Planck equatiorPS! there- 
fore 

(3) 

where 

mQe) = - 1^ (/ dpe-^^(''-'3^^^ . (4) 

The quantity F{Qe) is the equilibrium free energy for 
the crystal-order parameter, /3 = l/k^T^ and D = 
{{SQe)'^) /tqq is the diffusion constant projected along 
the Qa direction. The quantity {{SQq)^) « 0.01 is the 
mean-square fluctuation of Qe in the liquid basin for the 
216-molecule system considered in Fig. |2ja). The long- 
time limit is set by the diffusion constant, 

lim Pn,{Qe,t) xexp[-f3F{Qe)]. (5) 

Dt—¥oo 

For quantitative treatments of the ultimate equilibra- 
tion (i.e., of the final stages of crystal coarsening), Eq. [s] 
could be generalized to include Qg-dependence and mem- 
ory effects in D. Such generalization could account for 
the complexity of pathways by which multiple ordered 
domains reorganize and connect and would be expected 
to increase the timescales for equilibration. These refine- 
ments are unnecessary for the current analysis of early- 
stage coarsening, where Qe does not progress far from its 
values in the liquid. 

We have integrated Eq. |3] using a first-order finite dif- 
ference approach with a small enough discretization of 
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a) p = 2.0kbar b) p = 2.2kbar c) p = 2.4kbar 
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Qe Qe, Q& 



FIG. 3. Non-equilibrium pseudo free energy surfaces, Fnc{p,Q(i,t), at three different pressures, illustrating how artificial 
polyamorphism arises as a finite-time effect. All three surfaces are evaluated by propagating from an initial liquid distribution 
for a time t — Wtq^ . Computed as in Fig. [2j with the same notation as used in that figure. Contour lines are separated by 1 
ksT and statistical uncertainties are about 1 fcaT. 



Qq and time to ensure numerical stability.'^ Figm'e [2] (b) 
shows how the distribution evolves in time from an ini- 
tial Gaussian distribution centered in the liquid region of 
Qq. What is notable is that the relaxation to equilibrium 
takes orders magnitude longer than the basic timescale, 

With this time evolved probability distribution, we 
have used Eq. [ijto estimate a non-equilibrium joint free 
energy, 

Fnc{p,Q(>,t) = -kBT \n[P{p\Qe) Pnc{Q6,t)] . (6) 

This pseudo free energy function for two intermediate 
times is shown in Figs. [2] (c) and (d). At the first of 
these intermediate times, t = 10 tq^, i^nc(p, Qe,^) ex- 
hibits two minima at low values of Qq, and these minima 
are separated by a small barrier of a few k^T. The low 
density basin is centered at the mean density of the crys- 
tal, and the high density basin is centered at the mean 
density of the liquid. We use the reduced density vari- 
able, p = {p — pxtO/Ap, to emphasize these connections 
to the crystal and liquid basins. 

This behavior shown in Fig. l2] (c) is precisely the be- 
havior found in Refs. lUl El II of and [1^1- both the bi- 
stability and the length of time allowed for equilibration. 
Those workers find tq^ ~ 10® Monte Carlo sweeps, and 
they use 10^ sweeps to estimate free energies. The low- 
density liquid minimum eventually disappears, but the 
time scale for that to occur is orders of magnitude longer 
than considered in Refs. ITO l ll H ITS l andfTOl 

To further illustrate the connection between the non- 
equilibrium calculation shown here with the finite-time 
sampling results of Ref . UHl El |I3 and UHl Fig. |3] shows 
the effects of pressure variation on the pseudo free energy 
surfaces. Upon re-weighting to lower pressure, Fig. [3] 
(a), or to higher pressure. Fig. [s] (c), one of the disor- 
dered minima disappears. The specific dependence on 
re- weighting and the relative locations of the minima are 
also consistent with the results of Refs. jTUl El El and 
El We have thus reproduced the principal results of 



those papers by identifying the limited time over which 
the system was allowed to equilibrate. See for example. 
Fig. 2 of Ref. 15, as we have purposely used a similar 
color code in our Fig. [3] to emphasize the similarities of 
our finite-time results with the free energies reported in 
that paper. 

References El E] El and El employ slightly different 
variants of the ST2 model and slightly different temper- 
atures. In Appendix B, we consider another variant and 
another temperature to illustrate that the results of our 
analysis are generic. The analysis uses the Fokker-Planck 
equation in order to elucidate the generality of the phe- 
nomena. An explicit simulation calculation establishing 
the same conclusion is also provided in Appendix B. 



B. Phase flipping 

The time-dependent pseudo free energies illustrated 
above also shed light on previous reports of ph ase flip- 
ping between two seemingly distinct liquids .'^^'i^ In those 
reports, large transient density fluctuations occur inter- 
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t/lOOTp t/lOOTp f/lOOTp 

FIG. 4. Trajectories propagated with over-damped Langevin 
dynamics on the free energy surface pictured in Fig. [2|a). 
The trajectories are initiated in the liquid basin and run for 
insufficient times to pass to Q% values larger than 0.3. The 
trajectories thus illustrate early stages of coarsening in the 
ST2a model at T = 235 K and p = 2.2 kbar. 
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mittently between smaller amplitude motion while the 
system is globally liquid-like. This behavior is expected 
for trajectories of p when driven by the pseudo free en- 
ergy surface graphed in Fig. [s] (b). Indeed, Fig. |4] shows 
representative trajectories obtained by running over- 
damped Langevin dynamics^*^ on the free energy surface 
of Fig. [2ja) with the time-scale separation tq^ = 100 Tp. 
The trajectories were initiated in the liquid basin at tem- 
perature T = 235 K.l^ The trajectories look exactly like 
those presented as phase flipping in Ref s . [T5l and [TBI (The 
structural relaxation time in those studies is Tp > 1 ns.) 

Thus, so-called "phase flipping" is not a flipping be- 
tween distinct liquid phases. It emerges in a single liq- 
uid phase because there exists a separation of timescales 
between density fluctuations and long range order fluc- 
tuations. Large density fluctuations occur as the system 
is attempting to crystallize, accompanied by the forma- 
tion of sub-critical nuclei formation that shrink and cause 
further large density fluctuations. 

This behavior is transient, as the pseudo free energy is 
ever-changing with time. An analysis of time series would 
elucidate this nature of the phenomenon. Specifically, the 
mean transition time will show a dependence on trajec- 
tory length, reflecting the non-stationarity of the system. 
Consequently, the distribution of transition times will de- 
viate from simple Poisson statistics. Those offering phase 
flipping as evidence for two distinct liquids have not pro- 
vided a quantitative analysis supporting such statistics. 
Rather, studies illustrated in Ref. jlTj catalogue many 
independent trajectories some of which exhibit anoma- 
lously long quiescent periods between transitions lasting 
several times 100ns « 100 Tp. Data supplied for the ST2 
model in Ref. [T7] also illustrate that phase flipping oc- 
curs at pressures far below a proposed critical pressure, 
in apparent contradiction to the authors' claims that a 
critical point exists. The model behavior we discuss is a 
consequence of liquid-crystal coexistence, which by sym- 
metry must not have a critical point and would explain 
the insensitivity of this behavior to pressure. 

The only quantitive analysis provided in studies of 
phase flipping is the calculation of bi-modal density 
distributions, which are then fit to a universal Ising 
form, from which supposed critical parameters are 
extracted j^^'^^ ' The barostat used in those simulations is 
well known to not reproduce correct equilibrium density 
fluctuations,!^ making such a calculation difficult to in- 
terpret. 



III. CALCULATION OF REVERSIBLE FREE ENERGY 
SURFACES 

The previous section emphasizes the significance of a 
time-scale separation and the pitfalls for simulation that 
result from it. In this section, we detail procedures that 
can overcome the problems inherent in that time-scale 
separation, and we apply the procedures to the ST2 
model to evaluate the liquid-liquid hypothesis in that 



case. In light of recent speculations j ^^ l ^^ l ^^ * ^^ * we also 
analyze how these equilibrium surfaces are affected by 
changes in model parameters and boundary conditions. 



A. Methodology 

We have carried out Monte Carlo simulations with 
constant number of molecules, N, pressure, p, and 
temperature, T. This ensemble is appropriate for 
determining conditions of phase coexistences and is well 
suited for sampling dense liquid and crystalline states. 
Density, p = N/V, fluctuates with p and TV fixed because 
volume, V, fluctuates. The alternative Grand Canonical 
calculation for sampling global density fluctuation^^! 
can suffer from poor acceptance ratios (e.g., acceptance 
ratios are reported^" to be as low as 10~^ and can also 
be ill-defined in applications to crystals.!^ 

Trajectory algorithm designed for the task. The 

move set we employ for our Monte Carlo calculations 
is chosen to mitigate long correlation times expected 
for supercooled liquids and coarsening crystals. Specif- 
ically, to allow for collective reorganizations, we use a 
hybrid Monte Carlo algorithrrP^ that propagates an ini- 
tial configuration with Boltzmann distributed velocities 
under symplectic, norm preserving, molecular dynam- 
ics. For the models with internal degrees of freedom we 
use the SETTLE integrator while for single-site models 
we use a velocity Verlet integrator.'*^ The configuration 
is integrated for a time n St, where St is the integra- 
tion timestep. Each move is accepted with a Metropolis 
criterion,!^ so energy need not be conserved and conse- 
quently St need not be small. 

In practice, we generally range from St ^ 5 — 30 fs 
and the number of steps in a trial trajectory, n, to vary 
between 1-20 depending of the steepness of the free 
energy landscape. The choices are made systematically 
to minimize correlation times. Volume moves are used 
at a ratio of 2 hybrid Monte Carlo moves to 1 trial 
volume displacement. The relatively large value of St 
serves to swiftly propagate dynamics over long time 
scales. We have found that at supercooled conditions 
for N « 200, it significantly reduces correlation times 
for structural relaxation relative to energy conserving 
dynamics. For instance, in the case of the ST2 model 
discussed in detail below, the characteristic structural 
relaxation times under these moves are between 10^ 
to 10"^ Monte Carlo steps, depending on the specific 
value of density and temperature. In contr ast, single 
particle Monte Carlo moves reported previously ** * ** * * 
yield structural relaxation times that are between 10^ 
to 10^ Monte Carlo steps, and molecular dynamic^^^ 
yields structural relaxation times that are between 10^ 
to 10^ integration steps. Accounting for the factor 
n = 0(10), we see that our choice of hybrid Monte Carlo 
moves is computationally more efficient by 1-3 orders of 
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magnitude over single particle moves and by 2 orders of 
magnitude over molecular dynamics. This remarkable 
speed up must in part reflect the highly non-linear and 
correlated nature of dynamics at supercooled conditions. 

Controlled biasing methods. Fluctuations that re- 
sult in phase transformations are exponentially rare at 
conditions of coexistence or modest supercooling. Stan- 
dard methods of umbrella sampling get around the rare- 
event problem by adding biasing potentials, W{x^), to 
the Hamiltonian in order to enhance occurrences of oth- 
erwise improbable fluctuations. Re-weighting configura- 
tions correct for the biasing. Here, stands for a point 
in the configuration space for the system. The form of an 
added energy function can be chosen for convenience, but 
in order to guarantee the system will reach a stationary 
state it must be time-independent. 

Free energy methods that do not strictly adhere to this 
condition, such as meta-dynamics"*^ and Wang-Landau 
samplingjBSl converge only conditionally in the limit that 
the biasing degree of freedom is the slowest mode or when 
the change in the biasing term asymptotes to zero. This 
issue is particularly relevant for supercooled water be- 
cause pathways beyond the early stages of coarsening 
generally involve several slow variables in addition to the 
global crystal order parameter, Qq. An incorrect free en- 
ergy estimate will be obtained if one or more of those slow 
variables are not controlled in meta-dynamics or Wang- 
Landau algorithms. We do not exclude the possibility of 
inventing adaptive methods that could account for the 
physical behaviors of supercooled water and be more ef- 
ficient than the approach we employ, but such adaptive 
methods are not currently standard. 

The order parameters, p and Qq, are chosen to dis- 
tinguish phases of broken symmetries expected to result 
in water-like models at low temperatures. Our previous 
study demonstrated that Qg is a sufficiently sensitive or- 
der parameter to distinguish globally ordered from dis- 
ordered states accompanying a freezing transition.l^l This 
virtue noted, it must also be appreciated that Qq deviates 
from its disordered value only after a substantial amount 
of orientational order has developed in the system. This 
fact is demonstrated below and also in the Supplement 
to Paper 1^ 

The umbrella biasing potentials we employ are of the 
form 

l^(x^) = k [pix"") - p*] ' + At [QQix"") ~-Q*q]\ (7) 

where p{x^) and Qq{x^) are the density and crystal- 
symmetry order parameters, respectively for configura- 
tion x^ . The biasing potentials, with its force constants 
K and k, keep these order parameters close to their tar- 
get values p* and Qq . Each pair of target values defines a 
specific sub-ensemble or so-called "window" in configura- 
tion space. After collecting statistics in one window, the 
window is moved by changing the pair of target values, 
p* and Qg, whereupon statistics in the new window are 
collected. The procedure is carried out throughout the 



P-Qq plane, making sure that passage from one region to 
another is fully reversible, and that adjacent regions have 
sufficient overlap of statistics to enable further analysis. 

For this purpose, we find that k in the range of 500 
to 10,000 /cbT and k in the range of 1,000 to 2,000 k^iT 
cm^g^'^ is satisfactory. Statistics gathered in these bi- 
ased ensembles are then unweighted and the free energy 
differences between each ensemble are estimated using 
MBAR.IS21 

Stitching together data from different sub-ensembles 
yields the joint probability and thus free energy 
F{p,Qq;PtT), where p and T, respectively, denote the 
pressure and temperature at which the Monte Carlo tra- 
jectory is carried out. Assuming this function is deter- 
mined accurately over the relevant range of densities, free 
energies at other relevant pressures are determined from 
Eq. [2j i.e., by straightforward re-weighting.'^SUsni 

Care applied to ensure reversibility. In each win- 
dow, we initially equilibrated for up to 100 structural 
relaxation times as evaluated for the larger of either the 
initial liquid density or the equilibrium of bias window 
density. Production runs were between 50 and 1000 
structural relaxation times, depending upon the length 
of time required to obtain reliable statistics as judged 
from cumulative averages for Qq and potential energy. 
Here, structural relaxation time refers to the simulation 
time, t, required for mean square fluctuations in struc- 
ture factors to decay from their initial to 90% of their 
relaxed value. 

Three different techniques for generating initial con- 
ditions were used. Initial seeds were created by cooling 
an equilibrium liquid initially prepared at T = 330 K 
at a rate of 10 K/ns until it reached the target temper- 
ature. Seeds from this procedure were biased into dif- 
ferent windows in steps between adjacent windows, by 
gradually changing parameters of the biasing potential 
and with re-equilibration runs in between 
each step. At high Qq, the crystal that was sponta- 
neously formed using this procedure in all cases was a 
defected Ice Ic. For second generation seeds, we assumed 
that the spontaneously formed crystal was the relevant 
solid phase, so we prepared a perfect Ice Ic configura- 
tion, which was used to sample intermediate and high 
Qq states as well as bias them into low Qq regions to 
sample liquid states. Third generation seeds were ob- 
tained by melting an Ice Ic configuration and then using 
states along the melting trajectory to seed intermediate 
Qq windows. These configurations were subsequently bi- 
ased into the high and low Qq regions, again by gradually 
changing the parameters of W{x^), and again with new 
re-equilibration runs. In total, of the order of 10'^ inde- 
pendent biasing windows are used in the calculation of 
an individual free energy surface. Specifically, we gener- 
ated about 800 for each SW and mW free energy surface 
and about 2500 for each ST2 and TIP4P/2005 model free 
energy surface. 

One issue to keep in mind if one chooses to use parallel 
tempering and replica exchang^^U jg i]igx the shape of the 
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FIG. 5. Free energy surfaces, F{p, Qe',P, T), for three variants of the ST2 model at temperatures where others report evidence of 
liquid-liquid coexistence for the ST2 model. Possibilities of two-phase coexistence require changes in convexity, as re-weighting 
through Eq. [2]in that case can produce two basins of equal statistical weight. A coexistence pressure is then the value p + Ap 
at which there is equal statistical weight. For the three variants considered, the only changes in convexity are associated with 
coexistence between a liquid (low Qe) and a crystal (high Qe)- a) Free energy for the ST2a variant at T = 230 K and p = 2.2 
kbar, with N = 216. b) Free energy for the ST2b variant at T = 230 K and p = 2.2 kbar with = 216. c) Free energy for the 
ST2c variant at T = 235 K and p — 2.05 kbar with A'' — 216. See text for definitions of the different variants. Contour lines 
are separated by l.bksT and statistical errors over the surfaces average to less then 1 ksT. Quantitative features will change 
with system size. For example, as N grows, the mean value of Qg in the liquid basin will vanish as l/N^^"^ , while in the crystal 
basin it will remain finite. 



free energy surface changes with temperature. A liquid 
basin exists for T > Tg, but it does not exist for T < Tg. 
Moving between replicas that traverse this crossover will 
produce a transient shadow of the higher temperature liq- 
uid basin in the lower temperature replica. Appendix C 
illustrates how overlooking this behavior can lead to poor 
free energy estimates and false impressions of a second 
liquid basin. By using a control variable other than tem- 
perature, related methods might be designed to increase 
the relevant diversity of sampled configurations, but the 
underlying physics that makes a particular variable either 
applicable or inapplicable must always be considered. 

In principle, free energies can be correctly computed 
by any number of different methods, provided the pro- 
cedures are truly reversible. This property, reversibility, 
was explicitly checked in our calculations by constructing 
plots of all two-dimensional histograms and checking for 
hysteresis. Estimates of errors in free energy differences 
were made by computing overlaps and gradients of the 
distributions obtained by various routes. These steps, in- 
cluding bidirectional biasing to and from the crystalline 
phase and creating many independent realizations of ini- 
tial conditions, follow standard practices for computing 
free energies articulated in reviews such as Ref . UHl 

B. Results for different variants of the ST2 model 

Figure [5] shows free energy surfaces we have computed 
for three different versions of the ST2 model, variants 
that differ only in the manner by which long-ranged 
forces are computed. The phase behaviors in each case 
are similar, with one liquid basin and one crystal basin. 
Indeed, the existence of singularities in a partition sum is 
usually not sensitive to subtle changes in potential energy 



function. This is true because the existence of a phase 
transition is mostly dictated by dimensionality and the 
general form and symmetry of the potential energy func- 
tion.'^In contrast, the locations of the singularities (e.g., 
temperatures and pressures of coexistence) are often sen- 
sitive to subtle chan ges.^ Recent reports on variants of 
the ST2 modeP^I^ have hypothesized that differences 
in electrostatic boundary condition and non-electrostatic 
cutoff parameters can account for why our previously 
published resultJ^ found no liquid-liquid phase transition 
where others suggest it does exist. The results shown in 
Fig. [5] dismiss this hypothesis. 

The ST2a model. Panel (a) is for the model we used 
in Ref. [21 but at a temperature slightly lower than that 
considered in our earlier work. Our prior reported cal- 
culations were for T = 235 K, whereas Fig. [sf^a) is for 
T = 230 K. The results for T = 230 K differ very little 
from those at T = 235 K. The model employs a modifica- 
tion of the original ST2 potential for water. ^ The modifi- 
cation includes forces from the long ranged electrostatics 
that were neglected in the original model. The inclu- 
sion uses an Ewald summation with conducting boundary 
conditions. The non-electrostatic Lennard- Jones poten- 
tial is truncated and shifted at 7.5A. This model was re- 
ferred to as the mST2 model in Ref. |21 Here, we identify 
it as the ST2a model, and distinguish it from a modified 
ST2 model with insulating boundary conditions. 

The ST2b model. Panel (b) shows the free energy we 
have computed at the same temperature as considered in 
Panel (a), but now with insulating boundary conditions. 
In addition, a tail correction has been added to account 
for the portion of the Lennard- Jones potential neglected 
in truncation. Panel (b), therefore, shows our results 
for the equilibrated free energy surface of precisely the 
variant of the ST2 model used in Refs. [TUlandIT?! Here, 
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we identify it as the ST2b model. As in Panel (a), the 
surface exhibits a single crystal basin at large Qe and 
a single liquid basin at small Qq. The only significant 
differences between the two surfaces in Panels (a) and 

(b) are in the location of the liquid basin and the relative 
stability of the crystal. The density for the liquid basin in 
Panel (b) is higher than that in Panel (a), and the crystal 
stability in Panel (b) is reduced from that in Panel (a). 
Reasons for these differences will be discussed shortly. 

The ST2c model. Finally, Panel (c) shows our re- 
sults for the ST2 model using the variant described in 
Refs. [TT] and IT51 which we identify as the ST2c model. 
It uses a reaction-field treatment of electrostatic inter- 
actions and a Lennard- Jones tail-correction. This Panel 

(c) , plus Panels (a) and (b) together with the results of 
Paper I show that small changes in temperature and vari- 
ation in boundary condition have only marginal effects on 
the phase behavior of ST2 water. 

When effects of changing potential energy or temper- 
ature are marginal, the nature of those effects are easily 
interpreted and computed with the identitjPHl 



kBT ln{exp[-AH/kBT])p,Q, 



(8) 



Here, AF{p, Qq) is the change in free energy surface 
due to changing the temperature or Hamiltonian by the 
amount AH / k^T. The angle brackets refer to the ensem- 
ble average over configurations with the original temper- 
ature and Hamiltonian, and with the order parameters 
fixed at the values indicated by the subscripts. When 
ATi/k^T is large, or when its effects are large, calcu- 
lations with this formula will yield poor estimates of 
AF{p, Qg) because exponential averages are slowly con- 
verging.'^S] On the other hand, if AV. is small, or if its ef- 
fects are small, averages converge reasonably quickly. In 
that case, Eq. [8] becomes a computationally convenient 
estimator. We have checked explicitly that this measure 
of marginal behavior is satisfied with respect to the above 
cited variations for the ST2 model at T « 230 K and 
N « 200. For example, the typical size of AH /kBT for 
the comparison between Figs. [5] (a) and [s] (b) is 30% of 
the root-mean-square fluctuations in the net energy per 
kBT. 

The conducting boundary condition for Ewald sums, 
used for Fig. [sjja), strictly cancels electrostatic surface 
potentials in the energy function. These surface poten- 
tials arise from instantaneous polarization fluctuations 
and transient dipole moments of the total system. The 
conducting boundary condition is generally chosen for 
use with molecular dynamics simulations because the al- 
ternative, insulating boundary condition, typically intro- 
duces discontinuities in the potential energy whenever 
a molecule crosses the the periodic boundary.l^ For a 
dipole disordered system, like liquid water and ice Ih, 
this boundary condition should be irrelevant as its en- 
ergy will average to zero. 

Assuming the Ewald parameters are chosen such that 
energy is well converged in both cases, the pressure of 
the system with insulating boundary conditions is larger 



than that with conducting boundary conditions by an 
amount Apsmf . This pressure difference is given by^ 



Aps 



27r (2e - 1) ]VP 



(9) 



Here, e is the dielectric constant of the surrounding 
medium, M is the total system dipole, and V is the vol- 
ume. For a disordered system, the average M is zero in 
the thermodynamic limit, and Af ^ fluctuates with values 
that are a extensive in system size. Accordingly, Eq. 5 
shows that this strictly positive contribution to the pres- 
sure vanishes as ~ ^/V. More discussion on subtleties 
involved in the implementation of the Ewald summation 
is given in Appendix A. 

The non-electrostatic part of the ST2 potential is a 6- 
12 Lennard- Jones interaction. The simulations yielding 
Panel (a) truncate and shift this potential to zero beyond 
the oxygen-oxygen distance = 7.5 A. This truncation 
produces a net potential energy that is slightly higher 
than that with no truncation. An accurate correction to 
the potential energy that accounts for the neglected tail 
is the mean-fleld estimate —16Trej^,jcr^pN/3r^, where the 
Lennard- Jones energy and length parameters are eLj and 
a, respectively. Differentiation with respect to volume 
thus gives a tail correction for the pressure,'^' 



3r3 



Ap- 



'tail 



(10) 



Thus, while the two simulations yielding Panels (a) and 
(b) in Fig. [2] are both carried out at p = 2.2 kbar, the 
effective pressure of the former differs from the latter 
by the amount Ap ~ Apsu^f + ^Ptaib so that the mean 
density of the latter will be higher than that of the former 
by the amount 



Ap « Ap 



kBT (p)2 



(11) 



where the term in square brackets is the compressibility, 
{d{p) /dp)T- Evaluating Ap from the formulas above for 
the surface and tail corrections, and estimating the mean- 
square density fluctuations from the widths of the liquid 
basins in the free energy surfaces, we find Ap fa 0.1 g/cc, 
in harmony with the differences seen in Panels (a) and 
(b) of Fig. [5j The system is much less compressible in 
the crystal basin than in the liquid basin (i.e., the density 
fiuctuations are smaller in the crystal than in the liquid) , 
and as a result, the shift in position of that basin between 
Panels (a) and (b) is much less than that found for the 
liquid. 

The relative stability of the crystal basin in Panels (b) 
is notably less than that in Panel (a). This juxtaposition 
is another manifestation of the fact that with a given 
external pressure p, the effective pressure of the ST2b 
model is higher than that of the ST2a model. 

In the reaction-field treatment used to compute 
F{p^Q&) of Fig. [sjjc). Coulomb interactions are summed 
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directly up to a cutoff distance, Rc, and contributions 
from larger separations are approximated as those from 
an ideal polarizable continuum. With this approxima- 
tion, and assuming the medium has a large dielectric 
constant, the term 

N N 

- - 2 E E 'i?!^ [1 - (-^. - ^c)] , (12) 

i=l j = l c 

must be added to the potential energy evaluated with the 
truncated direct Coulomb sums. Here, fii is the dipole of 
molecule i, rij is the distance between molecules i and j, 
and the 0-function is unity for positive arguments and 
zero for negative arguments. This reaction-field approx- 
imation is reasonable for a homogeneous system,ESl and 
the asymptotic large dielectric assumption is isomorphic 
to the conducting boundary condition used with Ewald 
sums to construct Panel (a) . While, the potential energy 
computed with a reaction field method is not guaranteed 
to be the same as that computed by an Ewald summa- 
tion, judicious choice of cutoff in this instance results in 
reasonable agreement. Indeed, there is closer correspon- 
dence between Figs, [s] (a) and (c) than between Figs, [s] 
(a) and (b). 

The free energy surfaces at pressures other than those 
considered in Fig. [s] are easily produced by re- weighting, 
Eq. [2] Such re-weighted surfaces are shown in Appendix 
D. 



C. Global order contractiorJ^ 



p = 2.7 kbar 




FIG. 6. Free energy functions and mean order parameters 
computed for the ST2a model at T = 235 K, illustrating ar- 
tificial polyamorphism arising from incomplete knowledge of 
Qe-dependence in F{p, Qa). (a) The restricted contracted free 
energy function, F[p\ Q^^^), for A'' = 216 at several indicated 
choices of Q™''". The functions are computed from integrating 
the reversible free energy function in Fig. [2] (a) re- weighted to 
pressure p = 2.7 kbar. Error bars indicate one standard devi- 
ation. Higher or lower pressures shift the free energy to favor 
the liquid or crystal, respectively. See Eq. [2] (b) Re- weighted 
contracted free energy function, with Qg = 0.4, as if the 
re-weighting coincided with a Maxwell construction for two 
coexisting liquid phases, (c) The mean value of Qg as a func- 
tion of the maximum order-parameter value for T = 235 K 
and p = 2.7 kbar. 



Here, we explicitly demonstrate the importance of Qe 
for analyzing phase behavior of supercooled water by 
showing the effects of controlling the range of accessi- 
ble Qq fluctuations. This discussion follows closely that 
provided in Ref. [5^1 and it augments the results of Sec. 
II. 

In particular, we define a contraction of a constrained 
free energy, 

PF(p- gr'^) = - In dQe e'^^^^'^^)^ . (13) 

We compute these functions from our estimates of the un- 
constrained reversible free energy surface. Functions so 
obtained are shown in Fig. [6]ja). The unconstrained free 
energy from which they are derived is F{p, Qg) graphed 
in Fig. [2ja) and re-weighted to the pressure 2.7 kbar. It 
is the reversible free energy surface for the ST2a model 
at a pressure that puts the system close to coexistence 
between the liquid and the crystal at the temperature 
T — 235 K. At lower pressures, the liquid will be super- 
cooled; at higher pressures, the liquid will be stable with 
respect to the crystal. The surfaces for those different 
pressures are obtained from the pictured surface by ap- 
plying Eq. [2| 



An upper limit of Q™^^ ~ 0.13 encompasses the liquid 
basin. The contracted free energy in that case is uni- 
modal, and it does not exhibit statistically meaningful 
changes in convexity. That is to say, no re-weighting with 
Eq. [2]in that case will produce bi-modality. Thus, there 
is not a second liquid for all densities (and corresponding 
pressures) in the range considered. 

Notice, however, that this contracted free energy func- 
tion is skewed in a fashion where fluctuations towards 
low density are more probable than fluctuations towards 
higher density. In the limit of large system size, fluctu- 
ations within stable or metastable basins are Gaussian. 
The skewed behavior is therefore a flnite system-size ef- 
fect. Its physical origin can be resolved by increasing 
gmax p^^. Qmax _ 0.25, & shouldcr and change in con- 
vexity appears. For larger values of Q^^^, there is sys- 
tematic growth of the shoulder into a basin. This low 
density basin is the crystal. The mean value of Qg as 
a function of Q^^^, (<96)Qg'^>=: is shown in Fig. [6|c). 
This mean value remains at its liquid state value for 
gmax ^ 0.45. Thus, by sampling the full surface pictured 
in Fig.[2j[a), it is possible to identify the low-density basin 
as the crystal phase. 

When limiting the range of Qq to Q™'^^ < 0.45, the fea- 
tures shown in Fig. [6] could be easily misinterpreted as 



10 




FIG. 7. Reversible free energies and equation of state for 
liquid TIP4P/2005 at T = 185K and TV = 216 for pressures 
between 1 and 3 kbar. a) Free energy, F{p,Qq). Contour 
lines are separated by 1 kB.T and error estimates are less than 
1 fceT. b) Contracted free energy as a function of density 
when restricting Qe to the liquid basin, i.e., < 0.14. c) 
Liquid phase equation of state accessible from the free energy 
surfaces shown in Panels (a) and (b) 



indicative of liquid-liquid coexistence. Indeed, by apply- 
ing an external pressure to re-vifeight the curves in Fig. [6] 
(a), with Eq. [|as if F{p, Q™'''') 

was an equilibrium free 
energy for 0.13 < Qf"^"^ < 0.43, a bistable density dis- 
tribution can be obtained with a low mean value of Qq. 
This type of construction is illustrated in Panel (b) of 
Fig. [6[The bi-stability of precisely the sort reported in 
Refs. ITlIl [TTl [TSl and [T^ is thereby found. But contrary 
to the interpretation expressed in Ref. [I9l the behav- 
ior is not reflective of liquid-liquid transition. Rather, 
the unconstrained free energy shows that the appearance 
of convexity loss is associated with moving towards and 
then over the barrier separating liquid and crystal basins. 

As the geometry of the total F{p,Qq) changes with 
system size, N, the re-weighting (i.e., pseudo Maxwell 
construction) illustrated in Panel (b) and alluded to in 
Ref. [19] will depend upon system size in ways that are 
inconsistent with two-phase coexistence. The free energy 
barrier separating basins of truly coexisting phases scales 
as N"^/^. But the low density phase with Qg confined to 
low values cannot equilibrate and thus cannot coexist. 



IV. STUDIES OF ADDITIONAL COMPUTER 
SIMULATION MODELS 

The potential existence of a second liquid phase has 
been conjectured for a number of different models of wa- 
ter and other tetrahcdral liquids. In this section we sum- 
marize our work on applying robust free energy calcula- 
tions to a selection of these models. 



Larger phase space for the mW model 

In Paper 1,1^ we reported results for the mW model 
over a range of temperatures spanning 160 K to 320 K 
and pressures spanning 1 bar to 3 kbar. A corresponding 
states argumen1p2l indicates that this pressure range for 
mW water is a factor of 3 smaller than that accessed in 
experiments. As a consequence, the range of pressures 
relevant to the existence of a liquid-liquid transition in 
the mW model extends up to 10 kbar. We have thus 
expanded the domain of our free energy calculations and 
ruled out a liquid-liquid transition up to 10 kbar and 160 
K< T < 320 K. The free energy surfaces throughout are 
consistent with those shown in Paper I. 

Quantitative water model (TIP4P/2005) 

The recently parameterized TIP4P/2005 water 
modeP^ has had success in quantitatively reproducing 
many essential properties of water and ice, including the 
density versus temperature line at p = 1 bar. Previous 
studies have extrapolated equation of state data for this 
model to estimate the location of a putative liquid-liquid 
critical point and first order transition line.^^ Using 
free energy calculations we can test the validity of this 
extrapolation. 

Figure [7] shows F{p,Qe) for TIP4P/2005 computed at 
T = 185K and p = 1.8 kbar for N = 216 molecules. Ref- 
erence 56 estimates the location of the critical point to be 
T = 193 K and p — 1.35 kbar, as shown in Panels (b) and 
(c) of Fig. [7] The free energy in Fig. [t] is computed below 
this purported critical point temperature, yet the distri- 
bution is clearly monostable at low values of Qq. High 
values of Qg were sampled, but not shown for clarity. On 
re- weighting this free energy, Eq. [2j bi-stability does not 
appear throughout the pressure range, 1 kbar < p < 3 
kbar. We can conclude, therefore, that for this model 
of water, a second liquid does not exist at conditions 
studied. Here too, the putative liquid-liquid transition 
seems to be an artifact of finite-time sampling, reflecting 
a liquid-to-ice transition where coarsening is incomplete. 

General tetrahedral model (SW) 

We have also studied the behavior of the Stillinger- 
Weber (SW) model of silico n. Th is model has been the 
subject of numerous studie^^'Slli] f]^g^f^ have used equa- 
tion of state data to propose the existence of a second 
liquid phase. Using free energy methods, we can test 
this proposal by examining conditions below the puta- 
tive liquid-liquid critical temperature. 

Figures |8](a)-(c) shows F{p,Qq) calculation for a sys- 
tem of 512 particles at T = 1050 K and a range of pres- 
sures, < p < 10 kbar. Based on the phase diagram 
proposed in Ref. [131 also calculated with 512 particles 
and citing agreement with Ref. [SJ this temperature and 
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a) p = kbar b) p = 5 kbar c) p = 10 kbar 




FIG. 8. Reversible free energy, F{p, Qe), computed for the SW model for N — 512, T = 1050 K and three pressures. Contour 
lines are separated by 2 fee T and error estimates are less than 1 ksT. 



pressure range should traverse the first-order liquid-liquid 
phase boundary. Rather than finding two liquid basins 
upon decreasing pressure, we find the system crosses the 
line of liquid stability, Ts{p), for pressures lower than 
p = 1 kbar. 

Put into context with our discussion in Sec. II, the 
geometry of the free energy in Fig. [s] (a) illustrates how 
finite-time sampling can produce the illusion of liquid- 
liquid bistability. Specifically, the facile equilibration in 
density would result in an abrupt change between the liq- 
uid with density at p = 2.45 g/cc and an amorphous ma- 
terial with p = 2.35 g/cc, while over short timescales the 
slow diffusion in Qg would keep the value small. In fact 
this point was noted by the authors of Ref.[T4]when they 
state that out of 10 to 50 trajectories, "Non-crystallizing 
samples (an average of 5) were run for up to 10 relaxation 
times when possible." 



V. SUMMARY OF CALCULATIONS ON PUTATIVE 
LIQUID-LIQUID TRANSITIONS IN SUPERCOOLED 
TETRAHEDRAL FLUIDS 

Table |TT] summarizes all that is now known about the 
putative liquid-liquid transition in supercooled water and 
related systems. We see that the low-temperature be- 
haviors of several different models of water-like liquids 
are similar. The reversible phase behavior in all cases 
is that of one liquid, a liquid that can coexist with and 
transform to a lower density ice-like phase. In all cases, 
the time scales for density fluctuations in the supercooled 
liquid are several orders of magnitude shorter than those 
for long-ranged order fluctuations, and both time scales 
are strongly temperature dependent. 

These features lead to a rich non-equilibrium behav- 
ior, some of which we have illustrated here in our treat- 
ment of the early stages of coarsening, and previously 
in our treatment of confined supercooled water-f^^l This 
non-equilibrium behavior is clearly responsible for the 
numerous reports of polyamorphism in computer simula- 
tions of water, and it is likely important in processes that 



lead to the formation of glassy phases of water. While 
the former represent artifacts of finite-time sampling, the 
behaviors of glasses represent an important class of phe- 
nomena worthy of future study. Such far-from equilib- 
rium behavior, however, is beyond the scope of what we 
have considered here. 

The generic nature of what is found in so many mod- 
els makes it seem unlikely that plausible models of water 
will exhibit liquid-liquid coexistence and a second low- 
temperature critical point. It also suggests that any 
model exhibiting a crystal phase with lower density than 
the liquid phase will also exhibit transient behavior that 
looks like polyamorphism when viewed on the time scales 
no longer than those of early-stage coarsening. The rea- 
son why this non-equilibrium behavior is not observed 
in a recent study of the SPC/E model^^ is because that 
studj^ quenches from temperatures significantly higher 
than Tg for that model. 

The generic nature also implies that simplified mod- 
els, like mW water, can reliably serve as useful descrip- 
tors of water. Hypotheses on the nature of low tempera- 
ture water could often be studied in that way. Further, 
computational efforts that struggle with overcoming the 
separation of time scales inherent in low-temperature wa- 
ter could often be tested with such models. In our own 
work (see Appendix C below), it seems clear that er- 
rors in prior announcements of liquid-liquid transitions 
in supercooled water could have been detected by first 
examining the behavior of mW and SW models. 

In documenting the necessity of attending to time 
scales and relaxation, we have outlined robust methods 
by which equilibration and reversibility can be achieved. 
There is clearly need for independent assessment of this 
growing body of work, which we look forward to seeing 
in the future. 
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TABLE II. Summary of models and conditions considered in 
this and in previous studies 



Model 



AF'^T/K.p/kbar) FE'^ (r/K,p/kbar 



mW 
ST2b 

ST2c 
ST2a 
SW 

TIP4P/2005 



(242, 1.8f](240, 1.8 
(238, 1.9n(235, 2.0] 
(228, 2.2TT(224, 2.3 



(245, 1.8Q(240, 2.0fl 
(235, 2.2fl(230, 2.4f 



(1070, O.f] (950 7.5f] 
(920, 11.3|] 



(195, 1.45^ 



(160-300, 0-10.0|^ 
(230-240, 1.0-3.0 

(230-240, 1.0-3.o|j 
(230-240, 1.0-3.0|] 

(1050, 0-10.0 
(180-190, 1.0-2.6|] 



^ Conditions where artificial polyamorphism has been reported. 
^ Conditions where Qe-equilibrated free energy calculations rule 

out a liquid-liquid transition. 

Short Grand Canonical simulations from Ref. 1101 N 200. 

Qg-unequilibrated free energy calculations from Ref. 1151 

N 200. 

" Qe-unequilibrated free energy calculations from Ref. 1191 
N ^ 200 

^ Short molecular dynamics trajectories from Ref. 1141 N = 512. 
s Short molecular dynamics trajectories from Ref. 1561 N = 500. 

This paper and Paper 1, N = 216 to Af = 1000. 
' This paper, N = 216. 
J This paper and Paper I, N = 512. 
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Appendix A: Checks on coding and long-ranged force 
evaluations 



While the main text has postulated reasoned explana- 
tions for a discrepancy between the calculations of our 
pervious work and that of Ref. [TOl [11] [151 and [TH it does 
not address the possibility of an undetermined error in 
implementation. To address this we have worked with 
the authors of Ref. [15] in exchanging data and molecular 
configurations.'^ The results of these studies were unam- 
biguous: all parties were able to reproduce energy evalu- 
ations and standard mean and fluctuation quantities for 




100 200 300 400 500 600 700 800 

MC Sweeps/ lO'"* 

FIG. 9. Time series illustrating how some choices of param- 
eters for the Ewald sum in the ST2a model, which uses con- 
ducting boundary conditions, can lead to the formation of 
dipole ordered ice, while for other choices the liquid remains 
stable. See text. 



the ST2b model. Specific statistical properties computed 
were the average density and compressibility for a sys- 
tem of 200 molecules at two state points, T — 300 K, 
p = 1 bar and T = 235 K, p = 2.2 kbar. Agreement 
between calculations was found to be within statistical 
error.fS^ Moreover, tests were done swapping configura- 
tions of = 200 molecules and obtaining independent 
evaluation of the energy and Qe values. These tests 
demonstrated that both groups were evaluating the same 
energy function and the same order parameter. Thus, an 
implementation error does not seem to be the root cause 
of this discrepancy. 

As reported in Ref. T^, at specific regions of state 
space, simulations of a stable low temperature liquid are 
sensitive to details of the Ewald summation. Specifically, 
the authors of Ref.[Tn|found that with conducting bound- 
ary conditions and N « 200, liquid configurations spon- 
taneously evolve into a dipole ordered form of ice VII 
when pressure is elevated and T < 235 K. While the au- 
thors of that study report they were unable to find Ewald 
parameters under conducting boundary conditions that 
did not display this pathological behavior, other studies 
applying this boundary condition with other electrostatic 
estimators have not found this same difficulty. For exam- 
ple, Ref . [TH employs a reaction-field treatment of electro- 
static interactions with conducting boundary conditions. 

By juggling the method of evaluating Ewald sums, we 
too have been able to reproduce the pathological behav- 
ior reported in Ref. 1151 Two of our time series are shown 
in Fig. [9j Using hybrid MC dynamics for 216 molecules 
and conducting boundary conditions at T = 228 K and 
p — 2.2 kbar, we find that changes to the details of the 
Ewald parameters can result in either stable liquid be- 
havior (blue curve), or spontaneous dipole ordering (red 
curve). In the stable liquid case, the Ewald parameters 
are chosen to accurately estimate the energy and forces 
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p = 2.0kbar b) p = 2.2kbar c) p = 2.4kbar 




Qe Qe Qe 



FIG. 10. Nonequilibrium pseudo free energy surfaces, Fnc{p,Q6,t), at three different pressures, illustrating how artificial 
polyamorphism arises as a finite-time effect. All three surfaces are evaluated by propagating from an initial liquid distribution 
for a time t = Wtq^ . Computed for the ST2b model from the Fokker-Planck analysis with the underlying reversible free energy 
surface shown in Fig. [Sja. The appearance of artificial polyamorphism is like that found in Ref. |T5|for similar conditions. The 
time t = 10 TQg corresponds to the time used for averaging in Ref. 1151 Contour lines are separated by 1 ksT and statistical 
uncertainties are less than 1 ksT. The correct, reversible free energy surface for p = 2.2 kbar is shown in Fig. ISTb). 



of the long range part of the potential to 1 part in 10^ 
using a standard error estimator and a spherical wave- 
vector cutoff.^^ In order to maintain this level of accuracy 
with a changing box size the parameters are updated over 
the course of the trajectory. In the unstable trajectory 



ST2b, 4;„, = IOtq, 

T = 230 K, p = 2.2 kbar 




0.04 0.06 0.08 



Qg 

FIG. 11. Nonequilibrium pseudo free energy surfaces, 
Fna{p, Qsjt) illustrating how artificial polyamorphism arises 
as a finite-time effect. The correct reversible surface is shown 
in Fig. [5|b). The incorrect surface shown here is computed 
for the ST2b model from a non equilibrium distribution for 
Qe, where the distribution is obtained from simulations ini- 
tiated in the liquid basin and propagated for a time 10 tq^ . 
The surface shows the appearance of artificial polyamorphism 
like that found in Ref. [T^ for similar conditions. The simu- 
lation time tsim = lOTQg corresponds to the time used for 
averaging in Ref. 1151 Contour lines are separated by 1 ksT 
and statistical uncertainties are less than 1 ksT. 



a cubic wave-vector cutoff is used, and the parameters 
are held fixed which reduces the accuracy of the poten- 
tial estimate. There may be other ways to induce this 
pathological behavior, and these ways will depend upon 
the type of dynamics used. The principal point is that 
it is a finite-size effect that is sensitive to technical but 
unphysical details in implementing Ewald sums, and the 
effect can be avoided when sufficient care is applied to 
the algorithm. 



Appendix B: Early stage coarsening and artificial 
polyamorphism in the ST2b model 

The calculations carried out in Section II for the ST2a 
model have been also done for the ST2b model at T = 230 
K with N — 216 and pressures ranging from 2 kbar to 2.4 
kbar. These are the conditions and model considered in 
Ref.fTSl The results are shown in Fig. [To) The theoretical 
results agree with those found in Ref. 15, thus indicating 
that the free energies reported in that work suffer from 
finite-time effects. 

In addition to the Fokker-Plank analysis, we have com- 
puted PnciQeit) directly from a molecular simulation of 
the ST2b model, starting from an ensemble of liquid 
configurations and running for igim = lO^Qg. We have 
then convoluted that non-equilibrium distribution with 
the equilibrium P{p\Qq) obtained from the reversible free 
energy in Fig. [5|b) according to Eq. [I] The results are 
shown in Fig. |11| Again, behavior like that reported in 
Ref. Uni is obtained, thus indicating that the results of 
that paper suffer from finite-time effects. 
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Appendix C: Artifacts of un-equilibrated initial conditions 

A somewhat different but related source of systematic 
error would be found in simulations that do not com- 
pletely equilibrate from initial configurations taken from 
a higher temperature. For example, suppose one has es- 
tablished good equilibrium statistics for a system at tem- 
perature T H- AT, and then wishes to use configuration 
from that data to seed a free energy calculation at the 
temperature T. (One can consider parallel temperin^^H 
as one realization of this idea.) If at T 4- AT the system 
exists in a metastable liquid that becomes unstable at 
the lower temperature T, then configurations from the 
higher temperature will bias the statistics of the lower- 
temperature system for time scales short compared to av- 
erage time for Qg to equilibrate (i.e., to leave the liquid 
region). Therefore, if the runs performed at temperature 
T are too short, a remnant of the meta-stable liquid basin 
at temperature T -\- AT will remain in the estimate of the 
free energy at temperature T. Figure [12] illustrates this 
problem for the case of the SW model. 

The correct reversible free energy surface for the SW 
model at this condition is shown in Fig.jsjjc). In Fig. 12 
we show an incorrect surface, which is obtained by using 
equilibrated data at temperature T + AT — 1100 K as 
initial conditions, and then carrying out trajectories at 
T = 1050 K that run for only tsira = ^^tq^. Histograms 
produced by this procedure yield a pseudo free energy 
with a low-density liquid minimum at conditions where 
the actual liquid is thermodynamically unstable. Error 
estimates performed from the data collected in this way 
are very small, of the order of 1 k^T. These small error 
estimates reflect the highly correlated nature of the data. 
Changes of Qq are simply too slow to fully develop on the 
time scales probed. 



Appendix D: Pressure dependence of F{p,Qfi) for variants 
of the ST2 model 

References [15] and Uni suggest that the analysis of the 
ST2a model in Paper PI overlooks a liquid-liquid tran- 
sition because it examines a pressure that lies outside a 
hypothesized spinodal region. We examine the validity of 
this suggestion with Fig. [13] This figure shows the free 
energy surfaces at different pressures for the three vari- 
ants of the ST2 considered in the main text. The pres- 
sure variations are constructed by applying Eq. J2] to the 
surfaces graphed in Fig. [5] What is found in each case is 
that pressures higher than those considered in Fig.[5]shift 
stability towards the liquid and increases density of the 
liquid, and pressures lower than those considered in Fig.JS] 
shift stability towards the crystal and decreases density of 
the liquid. Variation of pressure over the ranges consid- 
ered by Ref s . fTSl and [T9l does not lead to liquid bi-stability 
in the reversible behavior of the ST2 model. There is one 
liquid phase, and no liquid-liquid coexistence or spinodal. 



SW, t,i„, = IOtq, 
T ^ 1050 K, p ^ 0.0 kbar 



2.50 




0.10 



Q 



6 



FIG. 12. Nonequilibrium pseudo free energy surface for the 
SW model obtained with initial conditions from T = 1100 K 
and short equilibration times. The correct reversible surface 
is shown in Fig. ]8]^a). Contour lines are separated by 3 fee T 
and statistical uncertainties are about 1 ksT. 
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